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Abstract 

Astrophysical sources of TeV gamma rays are usually established by Cherenkov 
telescope observations. These counting type instruments have a field of view of 
few degrees in diameter and record large numbers of particle air showers via their 
Cherenkov radiation in the atmosphere. The showers are either induced by gamma 
rays or diffuse cosmic ray background. The commonly used test statistic to eval- 
uate a possible gamma-ray excess is Li and Ma (1983), Eq. 17, which can be 
applied to independent on- and off-source observations, or scenarios that can be 
approximated as such. This formula however is unsuitable if the data are taken in 
so-called "wobble" mode (pointing to several offset positions around the source), 
if at the same time the acceptance shape is irregular or even depends on operating 
parameters such as the pointing direction or telescope multiplicity. To provide a 
robust test statistic in such cases, this paper explores a possible generalization of 
the likelihood ratio concept on which the formula of Li and Ma is based. In doing 
so, the multi-pointing nature of the data and the typically known instrument point 
spread function are fully exploited to derive a new, semi-numerical test statistic. 
Due to its flexibility and robustness against systematic uncertainties, it is not only 
useful for detection purposes, but also for skymapping and source shape fitting. 
Simplified Monte Carlo simulations are presented to verify the results, and several 
applications and further generalizations of the concept are discussed. 
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1. Introduction 

The field of very-high-energy (VHE, > 100 GeV) gamma-ray astronomy is 
currently in its third generation of instruments, with an advanced future project, 
CTA d, already being in its preparatory phase. The currently operated systems 
H.E.S.S. Q, MAGIC Q and VERITAS 0] have established more than 100 VHE 
sources in the skyS. The acceptence of these telescope systems covers few degrees 
in diameter, and declines smoothly with increasing distance from the pointing po- 
sition. Establishing a gamma-ray signal from an astrophysical source requires to 
significantly prove a gamma-ray excess over background events that typically re- 
main dominant even after selection cuts. This background is mostly composed of 
diffuse hadrons, part of which appears almost identical to electromagnetic show- 
ers and has to be considered to be irreducible |@]. Besides these "gamma-like 
hadrons", the irreducible background also contains smaller fractions of diffuse 
electrons and gamma rays. This irreducible background, and the statistical and 
systematic uncertainties that come with it, are one of the main limiting factors 
of TeV astronomy; usually, an observational campaign for a given source either 
reveals one source or none, and only in few cases, or if the effort of a large scan 
is undertaken, additional unexpected sources are detected. Therefore, a statistical 
source detection technique that is both sensitive to weak sources and stable against 
systematic effects is of crucial importance to the field. 

The standard test statistic to evaluate an excess of gamma rays from a given 
sky direction is Li and Ma [@], Eq. 17, which hereafter will be refered to as >S L m- 
It was established among several alternatives they evaluated in their paper, based 
on the fact that this likelihood ratio test statistic was the only one that yielded 
a satisfyingly Gaussian null-hypothesis distribution. The formula they presented 
was designed to compare the event numbers of an on- and off-source observation 
(N on , N ff), and allows for a scaling factor a between the effective observation 
times (/ on = a t ff) to account for unequal exposures. 

In modern observation practice, most Cherenkov telescope data are not taken 
in On/Off-mode (see Fig. [H left), because this strategy implies a lot of observation 
time dedicated to empty sky regions. Also, it is prone to systematic differences 
between the on- and off-data caused by instabilities in electronic or atmospheric 
operating conditions, especially if the off-data could not be scheduled contempo- 
raneously enough with the on-source observation. Therefore, usually the "wob- 
ble" technique B7J] is applied, in which the data are taken at two or more observa- 
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Figure 1: Observation schemes with asymmetric acceptance shapes (green areas). The sky posi- 
tion to be evaluated for a signal is marked with a yellow star, the corresponding off-regions with 
grey stars. Left: Original On/Off scheme. The off-data provides A^ff, the on-dataA^ n, and it needs 
one exposure ratio a = ? n/foff between those. Right: Wobble scheme with four observation posi- 
tions. The off data of each wobble set can be taken from all other wobble sets, resulting in four a 
parameters and their corresponding N on , N «- 



tion positions offset in different directions from the main target coordinates (see 
Fig. [0 right). In this way, each wobble set provides both on- and off-data at the 
same time. In its original idea, the exposure shape is considered to have some cir- 
cular symmetry, and for each wobble set the off-data can be taken from the same 
observation, at sky positions of similar distance to the telescope pointing position 
("reflected regions" [8]). In that case, the off-exposure ratio a is the same for all 
wobble data sets, and S LM can be applied after summing up the on- and off-events 
of all wobble sets. This constant a can also be achieved approximately if the 
off-data are taken from hadron-like background events ("template background", 
or a ring area around the source position ("ring background" | fl0(1 ). both 
of which require also symmetry assumptions or an efficiency correction through 
Monte Carlo simulations of the isotropic background. 
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In the general case, though, the acceptance symmetry might not hold and a 
Monte Carlo correction may involve too high uncertainties. This occurs for in- 
stance in very low-energy observations, where camera-based acceptance inhomo- 
geneities (dead pixels, trigger fluctuations) are both difficult to model in simula- 
tions and furthermore lead to features in the acceptance shape that can depend on 
the Alt/Az pointing direction or other operation parameters of the system. This 
is particulary troublesome in two-telescope systems like MAGIC yQ, where al- 
ready the basic geometrical overlap of the two fields of view implies an elon- 
gated, Alt/Az-dependent exposure shape. Under these conditions, the "reflected 
regions" approach does not hold, and the off-data for each wobble set have to 
be taken from the other wobble sets (see Fig. [TJ right). This approach can pro- 
vide a sensitive measurement, because with several wobble sets, the off regions 
are numerous and well-populated, but it results in a different a for each wobble 
set, which is not supported by Slm- On top of that, if this procedure has to be 
done separately for different types of data (be it for instance different Azimuth 
angles or telescope multiplicities), it leads to many more ar-parameters, and in 
general, some off-events may happen to be oversampled if they lie in more than 
one off-region, which is also not considered in Slm- As a consequence, while 
the background density can still be modeled under certain assumptions and some 
numerical effort [11], the test statistic Slm, if still applied in some way, becomes 
very approximative. This may be dealt with in practice by making additional high 
requirements to the signal-to-background ratio of a detection which however 
limits the effective sensitivity of the instrument. 



Therefore, unlike previous efforts BlOLIllh . this work will not pursue to extract 



the variables needed for Slm through the complex task of explicit background 
modeling. Instead, the likelihood ratio concept behind that test statistic will be 
generalized, and new formulae will be derived that can directly be applied to 
multi-wobble Cherenkov telescope data. To do so, no Monte Carlo simulations 
or exposure symmetry assumptions will be required. It will only be assumed that 
the telescope acceptance shape is the same for different wobble data sets if they 
are taken under similar operating conditions^. 

Besides that, this work will also address the disadvantage of 5 LM that it de- 
pends on the size of the signal region in which N on is calculated. This area is 
usually defined through an integration radius ("# 2 -cut", with 8 being the angular 
distance between reconstructed gamma direction and source position). Its op- 



Note that if this basic assumption does not hold, any other symmetry assumption also breaks. 
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timization either depends on the background density and an assumption of the 
source strength, or may involve several trials. In this paper, these assumptions are 
reduced by accomodating the knowledge of the point spread function (PSF) in the 
formulae, which makes them independent of the source strength or background 
density. 

Likelihood ratios are frequently used to convert a complex likelihood maxi- 
mization problem to a test statistic that follows a ^-distribution. This possibility 
was first proven in II 1211 . and was suggested for astronomical purposes in ft 1 311 . The 
technique is now widely used in counting type experiments, mainly in X-ray and 
gamma-ray astronomy, and can be applied both for detection and optimization 
purposes 111411 . It should be pointed out that criticism and potentially more accu- 
rate or more general alternatives to the likelihood ratio concept exist IU5U16U 17I1. 
but are not subject of this work. 

The structure of the paper is to first define the coordinate systems and naming 
conventions needed for the calculations. Then, a likelihood function is set up and 
maximized to gain all relevant parameters (Sec. [2]). Based on that, the test statistic 
is formulated in Sec. [3] and its intrinsic inclusion of S LM is demonstrated. In 
Sec. |H some further possible generalizations and applications of the formulae are 
discussed and a recipe for "Likelihood Ratio Skymapping" is suggested. Finally, 
some example toy simulations are shown in Sec.[5]to verify the method. 



2. Building the likelihood function 

In this section, a binned likelihood function is formulated which will be the 
basis of the likelihood ratio test statistic, but can also serve to fit the shape or 
positional parameters of a detected source. It is set up as a Poissonian probability 
function that evaluates the consistency of the different wobble subsets taken under 
a given operating condition with each other, allowing for a hypothetical source 
with a well-defined shape. It is a binned likelihood, but can be generalized to an 
unbinned likelihood easily (Sec. I4.5I ). 

In this paper, the term "operating condition" may in practice refer to a range in 
any quantity that influences the acceptance shape. This binning may for instance 
be done in the telescope pointing direction (Alt/Az) for two-telescope systems, or 
in the telescope multiplicity for a multi-telescope system like CTA. It may also 
be binned in atmospheric conditions, night sky background light level or discrete 
performance states of the camera. Of course, the formulae are equally valid if no 
such binning needs to be applied. 
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2.1. Namings and definitions 

The data are assumed to be taken as several wobble sets to = 1,2, ... ,W, 
centered at sky coordinates x obs >bJ ("pointing directions"), throughout operating 
conditions m = 1,2, ... ,M (see previous paragraph). Hence, the observation can 
be described as a set of WxM two-dimensional sky histograms. These histograms 
shall be set up in relative sky coordinates (x' = x - x b S ,a,) centered at the obser- 
vation direction x obs ^ of each wobble set (see grid on the "on-data" drawing of 
Fig. [T]). The individual relative sky bins are named i = 1,2,...,/, where i may 
be regarded as a representation of a two-dimensional bin (i x , i y ) without loss of 
generality. The number of gamma- like events in one such bin (to, m, i) is N^' m . 
In relative sky coordinates, the shape of the background event distribution is the 
same for all histograms that belong to a given operating condition, while a signal 
at a fixed absolute sky position will appear in different relative locations for dif- 
ferent wobble sets. In the following, the axes of the relative sky histograms (X' , 
Y' in Fig. []} are treated as synonyms for (relative) right ascension and declination, 
although in practice, one should replace those for axes that are truly rectangular 
throughout the field of view. 

The exposure of the data taken under a given operating condition m is dis- 
tributed among the wobble sets to in ratios of a™, such that Yiu a Z = 1- There are 
different technical solutions of how to evaluate these ratios, the simplest and most 
common being through the observation times 



'tot 

which is correct if the background rate is the same in all wobble sets of m. In the 
absence of a signal, the expectation value E(Nf ,m ) for a given bin (to, m, i) equals 
the background expectation nf'" 1 , which in turn is a well-defined fraction of the 
total background exposure taken under the operating condition m: 

£(Af m ) = nf m = cCn? . (2) 

If a signal is present, its shape can be modeled with a (normalized) kernel gf, 
which usually may be a function representing the gamma-ray PSF of the instru- 
ment, or a more extended shape for dedicated searches of extended sources. The 
gamma-ray PSF is usually a well-known performance characteristic that is deter- 
mined from simulations and/or observations of strong known point sources. It 
might, as any other data cut, slightly depend on the assumed spectral index of the 
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source. The kernel parameters gf depend on 10 because the relative source coordi- 
nates depend on the wobble offset. Also, the shape of the PSF (i.e. the resolution) 
might depend on the distance from the pointing direction. 

The kernel, which only characterizes the shape of the signal, is now multiplied 
by a scaling factor <p to constitute a signal term which is added to the expectation 
value as follows: 

E(N^\cf>) = aZnT(l + cfigf) (3) 

This way, the excess events implied by a given relative excess (p of a source is 
automatically proportional to the efficiency of the detector at the position i in 
the operating condition m, and to the exposure of the given wobble set co. This 
first-order approximation assumes that the background exposure nf' m , which is 
proportional to the efficiency to gamma-like hadrons, is also proportional to the 
gamma-ray efficiency (see Sec. l4.4l for a way to loosen this assumption). 

The expectation value of the relative excess cp from a gamma-ray source scales 
with its absolute flux. Therefore, estimators for <f> can be used for flux skymapping 
(after the efficiency correction outlined in Sec. 14.41) . Still, the meaning of (p is 
technically that of an excess (or deficit), and therefore it is not, like a flux, bound 
to be positive or zero. Also, its interpretation as a gamma-ray flux is not exclusive 
— a significantly non-zero excess can also be caused by other physical effects as 
for instance the moon shadow, or systematic detector artefacts not common to all 
wobble datasets. In any of these cases, the null hypothesis for the excess/deficit 
judgement hep = 0, which is a degenerate case of the signal hypothesis, and not at 
the edge of the parameter space. This is important for the following, because like 
this, the likelihood ratio test statistic can be expected to follow a ^-distribution. 
A remaining limitation of the cp parameter space is the fact that the expectation 
value for the event number, E(Nf' m \ cp), has to be positive, so (p has to be larger 
than-l/g^. 

While the normalization of the kernels gf has to be consistent among different 
a>, it neither has to be unity, nor does YjiSf have t0 be the same for different a>. 
For example, if additional off-observations are added to the analysis, they would 
have gf = for all L Also, Sec. 14.41 discusses a possible case where a varying 
normalization of gf is appropriate. 

For convenience in the following calculations, also the averaged kernel for a 
given operating condition is defined: 

sr = X<ur (4) 
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and the sum of events in a given bin (m, i): 

N m = Y j N?' m (5) 

at 

2.2. The likelihood function 

The likelihood function can be defined as the product of all Poissonian proba- 
bility functions of the bins in i, co and m: 

uk- i « •> nun w<, "i,^"" r " (« 



N ' \ 

m i iu i ' 



nn^^ wy (7) 



The latter step took advantage of the normalization of a™ within each m and the 
average kernel convention (Eq.|4]). 

For the maximization procedure, it is convenient to calculate the log-likelihood 
X = In L, which is 



(8) 



All terms that are independent of the free parameters n™, (p were absorbed into the 
constant K. 

2.3. Determination of the parameters 

The free parameters to be optimized in order to maximise the likelihood func- 
tion are I x M total exposures n" 1 of the operating conditions, and, if a signal at 
a given sky position is assumed, its corresponding (p. The following paragraphs 
show that this is possible analytically for the nf, and numerically for the relative 
excess parameter cp. 

2.3.1. Background density parameters nf 

To find the n'J 1 ^ that maximize the likelihood function for a given cp, one can 
calculate the first partial derivative of the log-likelihood function (Eq. [8]) for a 
given m! and i'\ 



dn m ' 

a> "t (9) 



8 



This expression equals zero if 



N m 



= TT7T=. (10) 



The second derivative of the log-likelihood function is always negative for this 
solution, so it is the maximum of the likelihood function for a given (f>. Note also 
that for the null hypothesis, = 0, the best approximator is, intuitively, 



<o = AT- (11) 



2.3.2. Relative excess parameter 

Inserting the optimized exposures (Eq. [TOl) into the log-likelihood function 
(Eq.[8]), the log-likelihood expression for a given (p parameter can be simplified: 



£(Nr\<f>) = f+lL 



-NT + ) Nf' n In AT 



(12) 



m i a) ' ' ' 

Again, ^' and represent terms that are independent of (p. 
The partial derivative w.r.t. is 

g£g|f> = y y y ^ sT^Z (M) 

As it turns out in simulations, finding the root of this expression, and verifying 
the positive second derivative, is numerically straight-forward. In all reasonable 
cases, it leads to one solution sup that maximises the likelihood function (see 
Sec. 14.2. 1 I for possible exceptions). 



3. The likelihood ratio test statistic 

3.1. The likelihood ratio 

The likelihood ratio is defined as the ratio of the maximal likelihood of the 
null hypothesis (0 = 0), and the global maximum of the likelihood, allowing for a 
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signal (0^0): 



A 



L{N?- m \n™ ) ,cp = Q) 

m i 

nnn(^. , 



n i,0 



,(1 + <t>gf)) 



(15) 



m i co 



The typical prescription to convert this to a test statistic |12D is TS = -2 In A. 
Since the null hypothesis has M x I parameters, and is a special case of the alter- 
native hypothesis, which has M x I + I parameters (because of 0), the difference 
in degrees of freedom is 1, and TS can be expected to follow distribution. 
As in the case of 5 LM , this only holds for high count numbers; in fact, since the 
underlying mathematics are the same, it can be assumed that the validity into the 
low-count regime is similar to that of 5 LM - 

Since VTS = -yjx^ is a half-normal distribution, 



-V^ X?? IX~I»(I^) (.6, 

is the Gaussian significance of the considered sky position to contain a gamma- 
ray excess (or deficit). This holds independently of whether the initial assumptions 
about the gamma efficiency and PSF are correct, since imprecise assumptions can- 
not contradict the null hypothesis and might therefore only make the test statistic 
less sensitive, but never wrong. Section 1531 shows a validation of this Gaussian 
nature in simulation. 

While the test statistic may be regarded semi-numerical due to the determi- 
nation of sup , it is unambiguous and can be determined at any desired precision. 
The trials usually needed to choose the optimal radius of the signal region are ob- 
solete, because the PSF of the instrument is incorporporated, which is typically 
well-known (see sec. |2.1D . 

For a given signal sup , the significance after Eq. [T6l depends on how much gf 
and g'" differ. Therefore, a higher number of wobble sets quite naturally leads to 
a higher significance without increasing the complexity of the analysis by manual 
selection of valid off-region(s). Furthermore, the formulae are invariant against 
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the number of different operating conditions M, as long as the coverage of the 
wobble sets is good enough to provide well-balanced exposure ratios a™ needed 
for the gf m to differ. 

3.2. The degenerate case of on-off-analysis 

In this section, instead of the general multi-wobble case with variable expo- 
sure shape and an arbitrary PSF, a case of only one operating condition (M = 1) 
is considered, with separate on- and off-target observations (W = 2), and a sim- 
ple step function kernel that represents a predefined on-target sky region. The 
observation times are t on = a t oS , so the above a^ relate to the a as follows: 

fli = ~~~~r (IV) 
a + 1 

a 2 = -J-r (18) 
a + 1 

Here and in the rest of this section, the index 1 refers to the on-observation and 2 
to the off-observation. The step function kernel is 1 within an arbitrarily shaped 
signal region (i e SR), and zero elsewhere. Consequently, the kernel constants are 



u _ 1 1 if to = 1 and i e SR 
gi ~ (0 else 

gm = fs+T if ^SR 
Si (0 else. 

Inserting these terms in Eq. [QJ it can be converted to 



(19) 
(20) 



d£(N?> m | 0) ZitSR N? =i "« ZfeSR Nr 2 

+ S TTTT^-TT (21) 



dtp (a +1)(0 + 1X^ + 1) (or + 1X05+T + 1) 

Won aN off 

(a+ 1)(0 + 1)(0^ + 1) (a + 1)(0^_ + 1) 

the root of which can analytically be calculated to be 



(22) 



= ■ (23) 

aN oS 



Putting this to Eq. [161 one finds 



A, on m (i^-^V) + N oS in ((1 + 



off, 



(24) 
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which is the well-known Eq. 17 of Li and Ma |@], making it a degenerate case of 
Eq. [16] of this work. The advantage of the latter is that it combines an arbitrary 
amount of differently populated wobble (or off-target) observations, and at the 
same time uses the actual PSF of the instrument instead of a discrete on-target 
area. This also avoids the need for an ambiguous choice of integration radius. 



4. Further generalizations and suggested applications 



4.1. Establishing several sources and their parameters 

Once a source is detected, the log-likelihood function (Eq. [13]) is fully ade- 
quate to estimate the parameters of the source, such as the position or extension. 
These source parameters can be regarded as parameters of the kernel gf (which 
the constant K" does not depend on). 

If one or several sources could be established with relative excesses <p n , and 
corresponding kernel parameters gf n , gf n , they can be inserted to Eq.[3]in order to 
scan the field for other sources, and cross-check the new null-hypothesis distribu- 
tion: 

E(Nf' n | </>) = aZ nf (1 + cf>gf + J] <t>* §V ( 25 ) 

n 

In this case, the log-likelihood derivative (Eq. [14]) turns to 

d£(Nr | 0) v. v. v. . _ gfd + Zn <Pn 8t) - «T(1 + 5L. *n O 



ZEE 



NT 



(26) 

and also the test statistic derived from the likelihood ratio is now different, because 
the null hypothesis already assumes the presence of sources: 



S = Vts = 



(\ + 



1 + 



sup g;" 



(27) 



Using these modified formulae is correct as long as the sources and their effective 
off-regions are spatially independent of each other. If not, then their averaged 
kernels gf n overlap, and the (p n cannot be optimized independently, in which case 
a multi-parameter maximization, or an iterative scheme has to be applied. 
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4.2. Likelihood Ratio Skymapping 

A common task with Cherenkov telescope data is to scan the whole field of 
view for unknown sources. Since the formulae provided here are perfectly viable 
to do this even in the general case of an unknown, operating-condition-dependent 
acceptance shape, I suggest the following "Likelihood Ratio Skymapping" proce- 
dure: 

1. The events are filled into histograms Nf' m , and the relative exposure ratios 
a™ are determined. 

CO 

2. A grid of absolute sky positions is defined. The bins are independent of 
those of the wobble histograms and should be large enough to avoid over- 
sampling (see Sec. 14.2.11) . 

3. For each grid point, the following procedure is applied: 

(a) The kernel constants gf are calculated for each wobble set to in relative 
sky coordinates. 

(b) The average kernels g™ are calculated, using the weights a™ (Eq.|4]). 

(c) The root of Eq. [14]is determined to obtain sup . 

(d) The significance is calculated after Eq. [161 using the sign of </> sup to 
distinguish positive from negative excess. 

4. If a source is found, it is modeled with a log-likelihood fit (Eq. [T3l) . and the 
whole procedure is repeated with the modified formulae of 14. ll until the sig- 
nificance distribution follows the null-hypothesis distribution of a Gaussian 
function. 

In this scheme, which follows a similar iterative strategy as the likelihood 
analysis of the Fermi Science Tools^l, the final result is a number of established 
sources, an empty significance skymap, and a significance distribution that meets 
the null hypothesis. There are several ways one might display these sources in a 
publication while avoiding the negative excesses elsewhere, but in this work, no 
fixed recommendation concerning this graphical issue will be pursued. 

4.2.1. Technical notes and caveats 

Calculating and interpreting a skymap requires some more considerations than 
the mere calculation of a significance value. Some of these concern features that 
can be regarded as caveats of the method introduced here, and some are features 
that equally occur in other VHE skymapping methods. 



3 http://fermi. gsfc.nasa.gov/ssc/data/analysis 
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1. Features common to most skymapping methods 

(a) In very poorly populated data sets, it can happen that the log-likelihood 
derivative (Eq. [14]) does not approach zero (in empty sky regions), 
or the likelihood increases towards negative <p until it meets its limit 
- i/£max- Both are caveats that equally apply to other skymapping 
methods, because the negative excess that can be considered is always 
limited by the fact that no negative event numbers can be measured. 
In Sec. |51 satisfying results are obtained by discarding skybins that 
belong to the former case (in areas without sensitivity), and using the 
limit of S when <p approaches -1/gmax- 

(b) If the grid of the skymap is small compared to the kernel/PSF of the in- 
strument, the single significance values S become correlated and their 
distribution cannot be easily tested to be compatible with a Gaussian 
or not. This oversampling problem is also common to most skymap- 
ping methods and needs to be addressed with an appropriate binning 
or trial strategy (see Sec. [5]), or an appropriate consideration of the 
correlations in the compatibility test. 

(c) In case of multiple sources (or, equivalently, very extended sources), 
signal photons of one source might appear in the effective off regions 
of another, weakening the significances of both. This issue can be 
addressed by observing in multiple wobble sets, and a sufficient sky 
coverage around a potentially extended source. 

2. Features specific to Likelihood Ratio Skymapping: 

(a) Since fewer assumptions go into the test statistic than in conventional 
background estimation techniques fl 1 ON - it is much less affected by sys- 
tematic uncertainties. However, it is also bound to be statistically 
somewhat less sensitive, at least in the two-wobble case where each 
of the two data sets only contributes one off region. If many wobble 
coordinates are used, every data set provides W - 1 off regions, and 
the reduced significance can be expected to be compensated rapidly. 

(b) As 5 LM , also the formulae presented here do not distinguish between 
excess and deficit of events. Therefore, the presence of a signal pro- 
duces a positive S peak at the source position, but also a negative S 
peak at the sky coordinates that the corresponding off-data are taken 
from. This effect cannot provoke false detections, and is compensated 
by the modified formulae in Sec. 14.11 Also, it looses importance with 
increasing number of wobble sets. 
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4.3. Excess events 

If a number of excess events has to be calculated, this can be done summing 
up the signal part of Eq. [3j 

A^ex = £ suP sr (28) 

m i oj 



In case the signal of other sources has to be considered, as discussed in Sec. 14. 1[ 
the corresponding expression is 

i ffi 

N ex = yy nt — . oo) 

4.4. Variable gamma/hadron acceptance ratio, flux skymapping and variable sources 

The relative excess parameter <p is based on the assumption that the gamma-ray 
acceptance is proportional to the acceptance for gamma-like hadrons. However, 
Monte Carlo simulations may indicate that this is not the case, and provide cor- 
rection factors 

= ^"(gammas) 
7i ef (hadrons)' ^ 

where the e are the efficiencies to gammas and hadrons. In the scheme presented 
here, these constants may simply be incorporated into the kernel constants: 

gt^Yl'gt (32) 

This way, or if y'" = 1 is a good approximation, the factor between hypothetical 
source flux and relative excess parameters sup will be constant throughout the 
field of view, and sup may be used to display the physically more relevant relative 
flux skymap (possibly omitting unphysical negative flux values for consistency). 
One has to bear in mind though that sup is a value that is relative to the background 
density, so it is only reliable if the background density is sufficiently high. In 
poorly populated skymaps, the number of excess events (Eq. [291) may be a more 
stable parameter. 

Another considerable case where the kernel normalization may be modified is 
when the signal is not constant in time, but variable, and an a-priori assumption 
of the light curve is available. 
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4.5. Unbinned analysis and "orbit" mode 

The presented formulae can be adopted for unbinned analysis if the limit of 
an infinitely fine sky-binning I is assumed. In this case, most bins are empty 
(Nf' m = 0) and the sums in Equations [13] JT4] and [16] turn into sums over the 
relative coordinates xj of the single events j = 1, 2, . . . , J (i.e. bins for which 
Nf' m =1). Instead of the previously discrete kernel constants, functions g^Cxj) 
have to be used. In this case, the log-likelihood function is 



£(x\\cp) = K" + £ln 



l+0^(xj) 
1 +0g m ;(xj)J 



(33) 



the derivative w.r.t (f> is 
and the test statistic is 



I 



^(xj)-g^(xp 



■^(l + 0^(xj))(l+0^(xp)' 



(34) 



2x^]ln 

7 



fl + 0, 



r^(x')^ 



sup 



1 + &upg m ;(xp 



(35) 



Although potentially more precise, this unbinned quantities may be computation- 
ally much more expensive, since the function g w (x') has to be calculated J x W 
times for each evaluation of X, whereas for a binned analysis, only I x W calcula- 
tions are needed, which are typically much fewer calls. 

A case where the unbinned approach might be computationally more effective 
than the binned approach is the application to "orbit mode" observations 1180, 
in which the telescope pointing is not discrete, but moves along a circle around 
the target coordinates. If one introduces the wobble angle ij that parametrizes 
the pointing direction along that circle, the above equations still hold if g Wj (xj) is 
replaced by a variable kernel g(xj, £), and g"^(xj) is the convolution of g(xj, £) with 
the distribution of £, within the operating condition m ; that the event j belongs to. 



5. Monte Carlo Simulations 

To verify the proposed Likelihood Ratio Skymapping scheme and its practical 
feasibility, two standard scenarios were simulated for which the method might be 
useful. For simplicity, again arbitrary X-Y-coordinates are used instead of right 
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ascension and declination. The histograms Nf' m are 3° x 3° and the kernel (PSF) 
is assumed to be a Gaussian function with cr 39% = 0.05°. The bin sizes are 1 cr 39% 
(although they may in principle be arbitrarily small, at expense of computational 
load). The sky grid is conservatively chosen to be very coarse ( V27r cr 39 % j^, to 
minimize the correlations between the significances. In practice, a finer grid may 
be chosen if the oversampling correlations are taken into account, which is not 
pursued here for clarity. The exposure ratios a™ are considered to be on-time 
ratios, and the total number of generated background events are 80000. For each 
of the simulations, the following plots are provided, both for in the presence and 
absence of a source in the field of view: 

• Measured distributions of events (Nf' m ). 

• Kernels gf and gf (Eq.|U), assuming the coordinates of the source. 

• Likelihood as a function of cf>, assuming again the coordinates of the source, 
to show how well sup is defined. 

• Resulting significance skymap and distribution. 

The root-finding was done with the bisection method, and no other approxi- 
mations were made, although the fact that most summands in Eq. [14] are almost 
zero may easily be exploited to build a more efficient algorithm. 

5.1. Case 1: W=2, M=l, asymmetric exposure shape, equal on-time, source at 
observation center 

This is, for an asymmetric exposure, a simple, ideal case. The source is at 
(0/0) in absolute sky coordinates, and the two wobble positions are offset by ±0.4° 
in X-direction. Figure [2] shows the results. No significant excess can be found and 
the null-hypothesis distribution is compatible with a Gaussian with a mean of zero 
and with a width of cr = 1. 

Despite the asymmetric exposure, S LM can be evaluated correctly, summing 
up the off-events for a given on-region in the same region of the opposite wobble 



4 The effective sky area covered by a two-dimensional Gaussian kernel is 2ncr 2 . In order for 
this to equal the area of a rectangular grid cell, and thus lead to the same number of trials as if the 
kernel had the shape of a grid cell, the grid constant has to be yfln cr 
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set. As mentioned i n | 3.21 the radius inside which events are integrated is basi- 
cally a free parameter so the significance is scanned in a range of radii between 
1-3 crp SF , to make sure that the comparison is not artificially biased in the favour 
of the new formulae presented here. In the absence of a signal (Fig.©, the signif- 
icance found at (0/0) is 0.4 cr with Eq.[l6]and S LM = 1.1 cr. 

Figure [3] is the same scenario, but inserting 300 excess events at the observa- 
tion center. The excess is detected at 6.6 <x with Eq. [16] and S LM = 6.3 cr. Since 
the source happens to be exactly on a grid point (at (0/0)), basically only one grid 
point shows the signal, proving that the correlations between the significances are 
negligible. In the significance distribution, this one value peaks out while not dis- 
turbing the general Gaussian appearance of the bulk of entries. As negative (p are 
allowed (in order to see the negative part of the distribution), the ambiguity of 
excess and lacking acceptance causes minor downward artefacts at (0.0/ + 0.8) 
in the skymap and at the negative end of the significance distribution. This is an 
expected feature of the algorithm (see Sec. 14.2. IT) . 

Finally, Fig. |4]is the same setup as before, but the detected source is included 
in the null hypothesis to validate the detection and look for possible other sources 
(see Sec. I4.ll ). The significance map is again empty, the downward artefacts are 
gone and the new null hypothesis is met. 

5.2. Case 2: W=4, M=2, pointing-dependent exposure shape, very different ex- 
posure times, one additional off-data set, extended source with large offset 

This is a very complex scenario that benchmarks the flexibility of the test 
statistic and the Likelihood Ratio Skymapping technique. The background events 
are distributed over 7 sub-histograms consisting of 3 wobble sets taken in two 
different operating conditions, and an additional set of off-data for one of the op- 
erating conditions. As can be seen in Fig.|5l the data can be combined without any 
problem, yielding an empty skymap and a Gaussian null-hypothesis distribution. 

Figure |6] is the same scenario, but with 1200 photons added to mimic an ex- 
tended source (cr src = 0.2°) at the large-offset sky coordinate (0.4/1.0). Even using 
the same PSF-Kernel as before, the source can be detected (5.5 cr). 

5.3. Gaussian property ofEq. \T6\and comparison to S'lm 

Figure [7] shows a simulation of significances calculated for the whole field 
of view of 20 signal-free skymaps simulated as in case 2 (Sec. 15.2b . The distri- 



5 There are recipes to determine which integration radius is appropriate, but a discussion or 
comparison of those is not within scope of this paper. 
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Figure 2: Case 1: W=2, M=l, asymmetric exposure shape, equal on-time, no source. In the upper 
left, the simulated measured distributions are shown, in coordinates relative to the observation 
center. Right of it are the kernels gf, assuming the source to be at the observation center, and 
below is the corresponding average kernel g™. The graph in the third row on the right is the log- 
likelihood as a function of <p, and in the bottom row, the skymap and significance distribution are 
shown (the dotted line is a Gaussian of /i = and cr = 1). The observation center is marked with 
a cross in the sky map, and green circles in the dala histograms. 
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Figure 4: Same as Fig. [3] but with the detected excess included as an established source in the 
null hypothesis. The significance map is again consistent with the null hypothesis. Consequently, 
the log-likelihood function peaks at a relative excess parameter ((f)) that is significantly larger than 
zero. 

bution contains also bins of comparably small exposure and bins for which the 
optimization of the ^-parameter reached its boundary condition. No treatment of 
the remaining correlation between the significances is undertaken. A Poissonian 
likelihood fit to the distribution yields Gaussian parameters of fx = -0.012 + 0.010 
and cr = 1.006 + 0.007. The^f 2 = 40.5 may be compared to the 37 non-zero bins 
of the distribution. This verifies both Eq. [16] as a Gaussian significance and the 
suggested recipes to deal with its features in the skymapping case. 

In the above case 1, the test statistic of Eq. [16] is somewhat higher than Slm, 
even though the latter involves an optimization of the applied integration radius 
(which is expected to lead to a selection bias). To study this possible trend sys- 
tematically, Figure [8] shows the distribution of many significances obtained for the 
same case simulated 1000 times. It shows that while Eq.[l6]does not always lead 
to a higher significance, there is a significant trend that in average it does. 

Repeating this study with a lower background level, however, the effect is 
much weaker, while with a higher background density, the effect is stronger. This 
is assumably because with a low background density, the optimized 6? 2 -cut will be 
large and include the whole signal, reducing the advantage of the known kernel 
from which Eq.[T6lbenefits. 

The conclusion can be drawn that the test statistic presented in this work per- 
forms slightly better in cases with low signal-to-background ratios, and similar in 
cases where the dominant statistical uncertainty is that of the signal itself. 
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Figure 5: Case 2: W=4, M=2, pointing-dependent exposure shape, very different exposure times, 
one additional off-data set, The example kernel plots refer to the sky position (0.4/1.0), marked 
with a cross in the significance skymap, where a source is injected in Fig. [6] 
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Figure 6: Same as Fig. [5] but with 1200 excess events introduced at the position marked with a 
green circle in the upper histograms. 
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Figure 7: Accumulated significance distribution of 20 simulated skymaps as in Sec. I5.2I De- 
spite the complex setup of several wobble positions and operating conditions, and the skymapping 
caveats outlined in Sec. I4.2.1I the significance distribution follows very well the shape of a Gaus- 
sian distribution (the solid line is the Gaussian fit to the data discussed in the text). 
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Figure 8: Comparison of significances obtained with Eq. [TBI and Slm (with optimized # 2 -cut). 
Despite the trials implied by the optimized 2 -cut, Eq. [16] yields a systematically higher signifi- 



cance. 
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6. Conclusions 

This paper derives a new generalized test statistic (Eq.[[6]) that can be used for 
significance calculation in VHE astronomy. The advantages over the existing test 
statistics are that it flexibly takes into account any number of data subsets from 
different wobble coordinates and operating conditions of the system, even if the 
acceptance shape is very irregular and different between these operating condi- 
tions. Also, it takes advantage of a known gamma-ray PSF while not requiring the 
optimization of an integration radius ("# 2 -cut"). 

The test statistic can be applied to any position in the field of view, so it is very 
suitable for skymapping purposes. The advantages of this approach is that the test 
statistic only makes minimal assumptions on the acceptance field of view and 
does not require any exposure symmetry or Monte Carlo simulations. It is hence 
unaffected by many systematic uncertainties. A "Likelihood Ratio Skymapping" 
procedure is suggested in Sec. 14.21 

The log-likelihood function (Eq. [T3l) can also be applied to fit the shape and 
position parameters of the source. The formulae are furthermore extendable to ac- 
comodate established sources in the field of view, a non-homogeneous PSF shape 
or gamma-to-hadron acceptance ratio, an unbinned analysis approach or the "or- 
bit" observation mode. If the background event density is sufficiently high, the 
relative excess parameter (p is well-suited to calculate a gamma-ray flux map of 
the field of view. 

In several simulated scenarios it is verified that the test statistic can not only 
handle the difficult situations it is designed for, but also seems to be systematically 
higher in the signal case, and therefore more sensitive, than the commonly used 
test statistic after Li and Ma [6], Eq. 17. 
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